start;

Lat = (-90:1:90)*glv.deg;
Lng = (-180:1:180)*glv.deg;

M = length(Lat);
N = length(Lng);
mag = zeros(M,N);

[X, Y] = meshgrid(Lng, Lat);
for ii = 1:M
    for jj = 1:N
        res = geomag([Lat(ii), Lng(jj), 0]);
        mag(ii,jj) = res(1);
    end 
end

h=surf(X./glv.deg, Y./glv.deg, mag);
set(h, 'EdgeAlpha', 0);
view([0,90]);
xlim([-180,180]);
ylim([-90,90]);
xylabelsEx('Lng / deg','Lat / deg');



